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Abstract 

We compute the B*Bn coupling for the static heavy-light meson using all-to-all 
propagators. It is shown that low-mode averaging with 100 low-lying eigenmodes indeed 
significantly improves the signal for the 2-point and 3-point functions for heavy-light 
meson. Our study suggests that the all-to-all propagator is a very efficient method for 
the high precision computation of the B*Bn coupling, especially in unquenched QCD, 
where the number of configurations is limited. 
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§1. Introduction 



Significant progress in B factory experiments has provided us with crucial information to 
determine the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements and to test the stan- 
dard model and the physics beyond. Among all the components, the element \V u b\ has 
attracted a great deal of attention with regard to the problem of testing the consistency of 
the unitarity triangle when combined with the angle sin(20i)[= sin(2/?)]. The consistency 
of the constraints derived from \V u b\ and sin(20i) is especially interesting because \V U b\ is 
determined by a tree level decay process, while the decay that determines sin (20!) has con- 
tributions from penguin amplitudes, which is sensitive to new physics. Although sin (20! ) is 
already known with a few percent accuracy, \V u b\ is known only within 10% - 20% from the 
inclusive and the exclusive semileptonic B decays. Moreover, the constraint on the unitarity 
triangle given by \V U b\ from the inclusive decay is only marginally consistent with that from 
sin (20i ). Therefore, it is important to reduce the error in the determination of \V u b\ from 
the exclusive decay, which can only be realized by improving the form factor calculation in 
lattice QCD. 

In general, extracting the form factors in lattice QCD is numerically more difficult than 
extracting decay constants or bag parameters, since it involves 3-point meson correlators with 
nonzero recoil momenta, whose treatment generally introduces non-negligible statistical and 
discretization errors. However, in order to determine the CKM matrix elements, such as \ V u b\ 
or \ V c b\, one does not need to know the form factor over a finite momenta, and in fact even the 
form factors at a single value of q 2 are sufficient. Thus, in principle, the form factors at zero 
momentum recoil can determine the CKM matrix elements, provided that the experimental 
data are statistically sufficient. This is indeed the case for the determination of \V c b\ from 
the B —> Dlv process. The situation for B — > nlu also seems promising, because there has 
been a rapid improvement in the quality of the experimental data that can be obtained. In 
fact, recently the BaBar collaboration observed the q 2 dependence of the B — > txIv decay 
quite precisely.^ 

It is known that symmetries often can be used to greatly simplify problems. In the 
case of B mesons, the approximate chiral and heavy quark symmetries restrict the hadronic 
amplitudes. In particular, the effective Lagrangian, which has both chiral and heavy quark 
symmetries, together with small symmetry breaking effects, allows us to obtain explicit 
relations among decay constants and form factors through systematic expansions in chiral 
perturbation theory and in terms of 1/M, where M is the heavy quark mass. An important 
example is provided by the form factors of B — > nlu semileptonic decay. By the chiral 
symmetry, these form factors can be expressed in terms of the heavy-light decay constants 



2 



and the coupling of the vector and pseudoscalar heavy-light mesons to the pion, Qb*b-k- In 
general the pionic coupling of heavy-light meson is defined as 



where cjq is the dimensionless coupling that appears in the heavy meson effective theory. 
Therefore, we can reduce the problem of computing the form factor to simpler problems of 
computing the decay constant and the coupling g b . 

The main sources of uncertainties in the types of computations considered here are (1) the 
quenching error, (2) the chiral extrapolation error for light quarks, and (3) the discretization 
and/or perturbative errors from heavy quarks. The first two of these are common to almost 
all quantities, and unquenched QCD simulations with quarks that are sufficiently light to 
solve these problems are being actively studied with various lattice actions. To reduce the 
third type of error, however, formulations of the lattice heavy quark which allow the non- 
perturbative renormalization and the continuum limit are necessary. Promising approaches 
have been proposed by the Alpha collaboration,® the Rome II group,®® and their joint 
collaboration.® In this approach, one uses the nonperturbative heavy quark effective the- 
ory (HQET) through order 1/M or relativistic QCD with a finite size scaling technique, or 
combinations of both. It has been found that in the first and third approaches, the results 
in the static limit play a crucial role in precise computation of the B meson. 

For the above reasons, there have been several calculations of the B*Btt coupling using 
HQET in both quenched and unquenched lattice QCD. In order to use B*Bir coupling for the 
precise determination of \V u b\, the sea quark effect from the unquenched calculation and the 
1/M correction are essential. The major drawback of HQET is that the static propagators are 
much noisier. This limits the accuracies and also makes it difficult to study 1/M corrections. 
In the quenched case, one can make use of a large number of gauge configurations to reduce 
the statistical error, while in the unquenched case, clever techniques to reduce the statistical 
error for a limited number of configurations are necessary. 

The Alpha collaboration® proposed using a smeared link for the HQET action, which 
can reduce the noise/signal ratio significantly. In particular, they showed that the HYP 
smearing® is the most efficient. 

Low-mode averaging is also known to improve the statistical accuracy. The TrinLat 
group® proposed a new method employing an all-to-all propagator by combining the low 
mode averaging and the noisy estimation. They showed that their all-to-all propagator 
significantly improves the heavy-light meson correlators, at least on coarse lattices. 

In this paper, we aim at a high precision computation of the coupling g^. With the goal 
of realizing high precision computations of the B*Bn coupling in dynamical simulations, we 
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carry out a systematic feasibility study on quenched gauge configurations, combining the two 
improved techniques of the HYP smeared link and the all-to-all propagators with low-mode 
averaging. We find that 100 low eigenmodes dominate the 2-point and 3-point static-light 
correlators after a few time slices. This implies that the low-mode averaging technique indeed 
does help to reduce the error significantly. Although the same error reduction is achieved by 
increasing the number of configurations in the quenched case, it has a big advantage in the 
unquenched case, where the number of configurations by the computational resources. 

This paper is organized as follows. In §2, we summarize the form factors of the B — > nlu 
decay and the implications of recent experiments with regard to the determination of \V u b\. 
In §3, we explain the soft pion relations of the form factors using the chiral perturbation 
theory for heavy-light mesons. In §4, we review the two techniques to improve the precision, 
i.e. those employing the lattice HQET action with the HYP smeared link and all-to-all 
propagators. In §5, we explain the methods to compute the coupling g^. Details of the 
simulation are given in §6. We present our numerical results in §7. Conclusions are given in 
§8. 

§2. Form factor 

The form factors for B — > txIv decay are defined as 
(ir(K)\u^b\B(p B )) = f + (q 2 ) 



mi, — mi 



(pb + KY - "" B 2 V 

+/V)^h^, (2-i) 

qz 

where the momenta p B and k n are, respectively those of the B and n mesons and q = p B — k n 
is the momentum transfer. The value of q 2 ranges from to g^ ax = (mg- m^) 2 ■ Ignoring 
the lepton mass, the branching fraction of the semileptonic decay B — > irlu is given in terms 
of f + (q 2 ) only by 
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(2-2) 



As seen from this equation, the hadron matrix elements (7r(k. K )\u / y IJ 'b\B(pB)) are given by a 
linear combination of the form factors f + (q 2 ) and f°(q 2 ), and hence one must solve the linear 
equation for two independent matrix elements in order to obtain the form factors. Because 
one of the matrix elements vanishes at zero recoil, the numerical accuracy of f + {q 2 ) near 
the exact zero recoil is extremely poor. For this reason, the extraction of the form factors is 
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based on the matrix elements with small but nonzero recoil momentum. This gives rise to 
somewhat large statistical and discretization errors. In fact, previous lattice calculations in 
quenched QCD EUS) anc [ j n unquenched QCD® - ^ 1 have yielded determinations of form 
factors with typically 20% accuracies, and it seems difficult to reduce the errors with the 
present technique. Fortunately, at zero recoil, one can use the soft pion relation to relate 
the form factor to simpler quantities, i.e. the decay constants and the B*Btt coupling, 

q 2 ^>m 2 m j n i — q /m B * 

Since the hyperfine splitting of B* and B is 45 MeV, the direct decay of B* — > Bit cannot 
be measured experimentally, but it can be computed using lattice QCD quite accurately. A 
quenched lattice computation of with HQET was first carry out by the UKQCD collab- 
oration, ® anc l subsequently Abada et al.® made a more precise calculation. Becirevic et 
al. also computed in rif = 2 unquenched QCD.® The currents value of are 

9oo = 0.51 ± 0.03 stat ± 0.11 sys for n f = 0, (2-4) 
9oo = 0.51 ± 0.10^ for n f = 2. (2-5) 

It should be noted that the quenched result is quite accurate and once one succeeds in 
obtaining the 1/M correction for the B meson either by interpolating the results in the 
static limit and the charm quark region or by computing the 1/M correction in HQET, the 
form factors at zero recoil can be predicted precisely. However, in unquenched QCD, the 
statistical error is one of the major sources of uncertainty. Given this situation, some clever 
techniques are needed to improve the accuracy with limited statistical samples. 

The experimental situation is quite interesting. The Belle and BaBar collaborations are 
measuring B — > nlu, and they have even measured the q 2 dependence. In the most recent 
results, BaBar measured the partial branching fractions of B — > txIv decay with twelve energy 
bins for q 2 of size 2 GeV 2 , where q is the momentum of the lepton pair: 

Ar(q 2 <q 2 <q 2 +1 ). (z = 1, ■ ■ • , 12) (2-6) 

We should be able to extrapolate the experimental data to g^ ax in the near future. Although 
the present experimental statistical error is still large, the super B factory is expected to 
produce 50-100 times more data, which would allow a precise determination of the differential 
decay rate in the soft pion limit, similarly to that for B — > Dlv decay. 

§3. Chiral perturbation theory for the heavy-light meson 

The low energy dynamics of the heavy-light meson and the light pseudoscalar meson 
system can be described systematically using the effective meson theory incorporating the 
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breaking terms of the chiral and the heavy quark symmetries. Defining the chiral field as 
£ = exp(iM//) with the light pseudoscalar meson field Ai and the heavy meson field H, 
which is defined in terms of the heavy-light pseudoscalar and vector fields B and B* as 



#=_(l + 7 %)[^ 75 + jB * 7 . 



(3-1) 



where the heavy meson field transforms linearly under the chiral and heavy quark symmetry 
transformations. The light-meson part is the ordinary chiral Lagrangian and the heavy-light 
meson part of the effective Lagrangian is 



-Tr [Hiv ■ DH] + gTr[HHA^ 5 ] +0(1/M), 



(3-2) 



where A^ = §(£^£ — 0^), and v is the velocity of the heavy meson. The pionic coupling 
of the heavy-light meson appears in the leading term of the Lagrangian. After including 

0(1/ M) correction, the effective coupling is g B * BlT = 9 + ^~~T7^~> where g\ and g 2 are the 

ivi 

coefficients in the next-to-leading term of the Lagrangian. The heavy-light V-A current is 
given by 



(V - Aft = -Tr[ 7 ,(l - l5 )HC] + 0(1/M), 



(3-3) 



where a = f B y/m B . Using this Lagrangian, it is straightforward to derive the soft pion 
relation systematically, including the corrections in the chiral perturbation theory, as well 
as 1/M corrections. At leading order in the chiral perturbation theory, and through order 
1/M in the heavy quark expansion]-^ we have 



9B*Btt{ 



(3-4) 
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where A Em g , — m g . 

The B*Btt coupling can be obtained from the form factor corresponding to the matrix 
element 



>B*(p B *,\)\A»\B(p B )) = 2m B *A Q (q 2 )—^q» + (m B *+m B )A 1 (q 2 
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(3-5) 



m B * + m B 

where p B and p B * are the momenta of B and B* and we have q = p B — p B *. The B*Bir 
coupling can be extracted from the residue of the pion pole with longitudinal polarization 

as 



9B* Bit 



2m B *A (0) 



[(m B * + m B )A 1 (0) + (m B * - m B )A 2 (0)), 



(3-6) 
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where we have used the kinematical constraint that the matrix element does not diverge at 
q 2 = 0. Thus the pionic coupling §q is given by 



m B * + m B ( . m B * - m B 
2^/m B m B » 2^/m B m B * 

In the static limit, Eq. (l3-7p simplifies to 

9q = 9oo = Ai(0). 



(3-7) 



(3-8) 



Now, the form factor Ai(0) in the static limit can be obtained from the matrix element at 
zero recoil as 

(B*\A t (0)\B) 



2m i 



Ai(0)e 



(3-9) 



which can be evaluated using lattice QCD calculations. 



§4. Lattice HQET and all-to-all propagators 



The lattice HQET action in the static limit is defined as 



l+7o 



h(x) - Ul(x - l)h{x - 4) 



(4-i; 



where h(x) is the heavy quark field. The heavy quark propagator is obtained by solving 
the evolution equation derived from the action with negligible numerical cost. Despite the 
advantages of being numerically economical and allowing better control over the systematic 
errors than other lattice heavy quark formulations, the lattice calculations of heavy-light 
systems with HQET suffer from large statistical noise. The reason for this noise can be 
understood as follows. First, consider the 2-point function of a heavy-light meson. The 
signal C*2(t) and the noise Z\C|(t) behave as 



C 2 (t) = ([Qq}(t)[qQ}(0))~exp(-E Qq t), 
AClit) = (\[Qq](t)[qQ](0)\ 2 ) - exp [-{E Qq 
Thus, the noise-to-signal ratio is 



E qq )t] . 



noise 
signal 



exp 



Enn + E, 



Qq 



(4-2) 
(4-3) 

(4-4) 



In the static limit, the corrections to the self-energy and binding energy of the heavy-heavy 
system exactly vanish, while in the heavy-light system, the power divergent self-energy cor- 
rection gives a significant contribution to the energy. Therefore, the noise-to-signal ratio 
is 



noise 
signal 



exp 



(5m Q + + ^)t 



(4-5) 
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where the mass shift has the regularization dependent power divergent term Stuq ~ c/a, 
with c being the regularization dependent constant. This implies that the noise problem 
becomes increasingly serious as the continuum limit is approached. 

In order to reduce the noise, the Alpha collaboration proposed a new HQET action,^ in 
which the link variables are smeared to suppress the power divergence. They studied the 
static heavy-light meson in the cases that the links are smeared by APE smearing, HYP 
smearing^ 1 and a one-link integral and found that the smeared link significantly reduces the 
noise. This result is indeed consistent with the observation that the power divergent mass 
shifts at one loop are reduced by smearing. For a f3 — 6.0 lattice with 5000 configurations, the 
noise-to-signal for the time extent of 1-2 fm remains at the 0.5 - 1% level. This means that 
using smeared links and data with high statistics allows precise computation of the static- 
light meson. Indeed, Abada et slW^ used the HYP smeared link action and obtained with 
160 configuration to 3-5% accuracy, while for the unquenched case with 50 configurations, 
the statistical error of g^ was on the order of 10%. 

Another useful technique to reduce the statistical error is to employ the all-to-all propa- 
gators using low-mode averaging. Low-mode averaging was introduced by DeGrand et alP^ 1 
for the overlap fermion and has been used extensively in the e-regime along with baryon 
propagators.^ The TrinLatP collaboration developed a comprehensive method using all-to- 
all propagators by combining the low-mode averaging and noisy estimator with time, color, 
and spin dilutions. Their method is as follows. First, define the lattice Dirac Hamiltonian 
Q as 

Q = l5 D, (4-6) 

where D is the lattice Dirac operator. Then, the quark propagator S q (x,y) can written in 
terms of the inverse of the Hamiltonian Q = Q' 1 as 

S q (x,y) = Q(x,y)-f 5 . (4-7) 

The hermitian operator Q can be decomposed into low-mode and high-mode parts by making 
use of the spectral decomposition as 

Q = Qo + Qi, (4-8) 

Qo = J2 X ^ (l) ® v(m ' Qi = Yl ® vim i (4-9) 

1=1 j = AT ov + l 

where Aj is the eigenvalue associated with the eigenvector with the index labeling eigen- 
modes in the order of increasing absolute value of the eigenvalue. Correspondingly, the 
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propagator Q is also decomposed as 

Q = Qo + Qi, (4-10) 
Q = ]T-^W 4)t , Qi= V —v^<B>v^ = QP!, (4-11) 

t=l j = JVcv + l J 

where 

P t = 1 - P = 1 - W W ® yWt (442) 

i=l 

is the projection operator that projects into the space consisting of only the modes above 
the iV ev -th lowest eigenmodes. Determining the low-lying eigenmodes, one can compute Qo 
directly and evaluate the remaining high-mode part, Qx, with a noisy estimator. 
Generating a random noise vector rj which satisfies the property 

((v(x)lv(v)f)) = ^jJ-E^C^wO/)? = M-6<U (4-13) 

and introducing with 

4> [r] (x) = J2(Q p i)( x >yHr](y), (4-i4) 

we can obtain the all-to-all propagator as 

(<Vfr](*)«® ^.](!/)&>> = (4-15) 

where we have used the property in Eq. (14- 13[) . The TrinLat collaboration® showed that the 
dilution is efficient for decreasing the statistical error of the correlator evaluated with noisy 
estimator. The dilution consists of the decomposition of the noise vector t] into each time, 
spin and color component as 

%](^,t)« = Z)^w( f '*)' t 4 ' 16 ) 

3 

where j is the index for the dilution labeling the time, spin, and color sources, i.e. j = 
(t Q ,a ,a Q ), and r¥){x,t) is defined by 

V${x,t)Z = V[r]{x) a J t>to S a>ao S am . (4-17) 
After applying the dilution, the all-to-all propagator reads 

N r 



{QPi)(x,y) = lim (4-18) 



N T —*oc 

r j 
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Table I. The quantum numbers of the mesons and the corresponding gamma matrices. 
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The TrinLat group showed that their all-to-all propagator is indeed useful for obtaining 
a good signal for the heavy-light system on coarse lattices. Although their proposal seems 
promising, because the power divergence becomes increasingly severe as the lattice becomes 
finer, the magnitude of the quantitative improvement is realized for finer lattices remains to 
be investigated. 

§5. Calculational methods 

The matrix element (B*\A fl \B) at zero recoil can be obtained from the ratio of the 3-point 
and 2-point correlation functions R(t) as 

<B ' ( y |B(0)) = lim H(t,t A ), (5-1) 
Zm B t,t A ^oo 

where 

m tA) ~ (o B *(t + t A )o B m ' (5 ' 2) 

with the operators O b and O b * having the quantum numbers of the B and B* mesons, 
respectively. We apply the smearing technique to enhance the ground state contributions to 
the correlators as 

B (t,x) = ^<j>{r)q{t,x + f) 7 5 /i(t,f), (5-3) 

r 

B *(t,x) = 5>(f)g(t,x + r) 1 i h(t,x), (5-4) 

r 

where (f>(x) is the smearing function. To apply the all-to-all propagator technique introduced 
in the last section, we decompose the 2-point and 3-point correlation functions into the high- 
mode and low-mode parts as described below. 
The heavy-light 2-point correlator is written 

C ^ = 7 E E ( Tr [ s *(y> °' f ' + f > t ~>y + ™> °) h^(w)] ) , (5-5) 

x,y f,w 

where Sb is the heavy quark propagator, and r and /g are given in Table [Q C(t) is 
decomposed into two parts: 

C(t) = (C (t)+C 1 (t)), (5-6) 



10 





t + t 



to 



t + t 



to 



Fig. 1. The diagrams for the 2-point correlation function. The left and right panels show the 
low-mode and high-mode contributions, respectively. 
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(5-7) 



-| ' 'CV -I 

C ^ = v E Y t E 5 Tr i 75 ^ (f ' t] ° m ~ f)75r 

=1 x,y r,w 

x [v {i) (x + f,t) ® v {i) \y + uj,0)] f B (f)(f)(j)(w) 

c ^ = vw r E E E Tr {^te *; o)5(y - f ) 75J r 

x +r,t) ® 77^(^ + ^,0)] f 5 0(f)0(u?)} • (5-8) 

Figure [I] schematically depicts the contributions from the low-modes and high-modes. 

While the above functions Co(t) and C\(t) are represented with a fixed-source time slice, 
t — 0, when the all-to-all propagator technique is applied, all the translationally equivalent 
correlators are averaged over. We note that one does not necessarily average over all the 
source time slices, because after the time dilution is applied, the source time slice can be 
chosen arbitrarily. The number of source time slices adopted should be chosen appropriately 
for each Cq and C\ by considering the statistical errors and numerical cost. In the case of 
the 2-point correlation function, we average over all the available source time slices for both 
the low-mode and high-mode parts. 

The heavy-light 3-point correlator is given by 

C(t B *,t A ,t B ) = ^EE( Tr S f>(y^B] x,t B *)r B *Q(x + f,t B *; z,t A ) 

x,y,z r,w 

xr Aa Q(z, t A ; y + w, t B )f B (f)(r)(f)(w) 

where .Tr* 
parts: 



(5-9) 

la, r Aa = 7 Q , and F B = — 1 (a = 1,2,3). They are decomposed into four 



C(t B *,t A ,t B ) — \Cu(t B *,t A ,t B ) + Cih(t B *,t A ,t B ) 

+Chh(t B * ,t A ,t B ) + Chi(t B * ,t A ,t B ) 



(5-10) 
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Fig. 2. The diagrams for the 3-point correlation function. The top left and right panels display the 
lower-lower and lower-higher parts and bottom left and right panels display the higher-lower 
and higher- higher parts. 

Here, Cu, Cih, Chh, and Chi denote the low- low, low-high, high-high, and high-low mode 
correlators, respectively, which are schematically depicted in Fig. [2 Explicit expressions of 
the 3-point functions are given in the Appendix. 

As noted for the 2-point correlation functions, translationally equivalent correlators can 
be averaged over an appropriate number of source time slices. For the low-low part, we 
averaged over all the source time slices, since the numerical cost is small. For the other three 
parts, we explore the best solution in the following section. 

56. Simulation details 



Our simulations were carried out on a quenched 16 3 x 48 lattice with the standard pla- 
quette gauge action (3 = 6/g 2 = 6.0. We considered 32 gauge configurations generated by 
the pseudo-heat-bath update algorithm, each separated by 1000 Monte Carlo sweeps, after 
a period of 1,000 sweeps allowed for thermalization. 

In the HQET, we used the static action with HYP smeared links and the parameter 
value (0(1,0:2,0:3) = (0.75,0.6,0.3). We used the C(a)-improved Wilson fermion for the 
light quark with the nonperturbatively determined clover coefficient csw = 1.769® We 
considered three hopping parameters, k = 0.1335, 0.1340, and 0.1342. The B and B* 
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meson operators at both the source and sink points were smeared with the smearing function 
4>(t, r) oc exp(— 0.45|f|). This smearing function was obtained from the wavefunction given in 
Ref.24) by properly scaling the parameters to our lattice spacing. The low-lying eigenmodes 
of Q were obtained with the implicitly restarted Lanczos algorithm,^ 1 accelerated by spectral 
transformation with the Chebyshev polynomial of degree 80P The low -mode parts of the 
correlation functions can be computed with these eigenvectors. For the high-mode parts of 
the correlator, we obtained the quark propagator with a source vector which is given by a 
normalized complex Z2 random noise and then projected into the space perpendicular to the 
space spanned by the low modes. Then, the time, color, and spinor dilution was applied, and 
all color and spinor components were summed over, while the number of the time dilution, 
N to , may be kept less than the total time N t = 48. The Dirac operator was then inverted 
using the BiCGStab algorithm with the stopping condition 10~ 10 . The number of eigenmodes 
for the low-mode averaging was chosen to be N cv = 100 and the number of the random noise 
and the time dilution was taken to be N r = 1 and N to = 4 for our final productive run. This 
choice was made so as to realize good statistical accuracy of the correlators while keeping 
the numerical cost at a reasonable level, based on our exploratory study of the 2-point and 
3-point functions by changing the parameter choices. This is explained below. 

6.1. Exploratory study of the parameter choice for all-to-all propagators 

In this subsection, we report the results of exploratory studies to determine the optimal 
parameter value for the all-to-all propagators in order to realize maximum numerical accuracy 
of the correlators with modest numerical cost. We found that obtaining 100 low eigenmodes is 
feasible, because it costs about two hours per configuration of our computational resources; 
this cost is not large compared to that of the other parts of the calculation. For the correlator, 
the numerical cost for the all-to-all propagators for the low-mode part grows linearly in 
N cv for 2-point correlators, since there is only one light quark propagator in the diagram, 
while for 3-point correlators it grows quadratically in N cv , because there are two light quark 
propagators in the diagram. Then, the numerical cost for the high-mode part grows linearly 
both in N to and N r , because we apply the source method to compute the high-high part of 
the 3-point correlation functions. 

In the exploratory analysis described below we studied the parameter dependence using 
the data with 32 configurations for both 2-point and 3-point correlators. All the studies 
of the 2-point correlator were carried out with N r = 1. Except for the study of the N to 

*' The computational costs cited in this section, are the CPU times using one CPU of our vector su- 
percomputer. However, the CPU time depends strongly on the architecture, and therefore a comparison of 
the cost does not have absolute meaning. Instead, if should be understood an only one particular reference 
value for guidance, among many others. 
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Fig. 3. Plots of the effective mass of the 2-point correlators for k = 0.1335, with low mode averaging 
and Nt = 4 and N r = 1 for high modes. The top, middle and bottom panels correspond to 
N cv = 0, 20, and 100, respectively. 

dependence, the 2-point functions were basically studied with Nt = 4, while in the final 
production run, we used N to = 48. For the 3-point correlators, basic studies were carried out 
for (N r , N to ,tA) = (1,4,8) except for the study of the (N r ,N to ) and tA dependences, while 
for the production run, the same parameter values were adopted. 

6.1.1. iV ev dependence of the 2-point correlator 

We first study the 7Ve V dependence of the 2-point correlators. The low-mode averaging 
gives the exact all-to-all propagator for the low mode part, and thus it improves the statistical 
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Fig. 4. Low-mode and high- mode contributions to the 2-point correlators for k = 0.1335 with 
low-mode averaging and Nf = 4 and N r = 1 for the high modes. The top, middle and bottom 
panels correspond to N cv = 0,20, and 100, respectively. 



accuracy of the low-mode contribution significantly. Therefore if the low-mode contribution 
dominates the 2-point function for the time range where the effective mass exhibits a plateau, 
one can expect that better accuracy can be obtained. Figure [3] plots the effective mass of 
the 2-point correlators for three choices of the number of low eigenmodes, N cv = 0, 20 and 
100, with N r = 1 for k = 0.1335. The statistical fluctuations apparently decrease as iV ev 
increases, and for iV cv = 100 we observe a clear plateau for t > 7. The situation is almost 
the same for k = 0.1342. 
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The behavior described above can be understood as the low-mode dominance in the 
correlator for large t. Because we used the low-mode averaging, the low-mode part has 
the best accuracy. For the high- mode part, we only averaged over N to = 4 timeslices for 
source points in this particular analysis, and for this reason, maximal statistical accuracy 
was not obtained. Also, the high-mode part suffers from additional statistical error due to 
the random source for the noisy estimator. These makes the high-mode part much noisier 
than the low-mode part. Figure. H] shows the low-mode and high-mode contributions to the 
2-point correlators for different choices of the number of low eigenmodes, N ev . We find that 
for N cv = 100, the low-mode contribution dominates the 2-point correlator for t > 6, where 
the plateau has almost been reached, while for 7V ev = 20, low-mode dominance is obtained 
only for t > 15. 

Although the 100 eigenmode calculation requiring 2 hours is rather time consuming, it 
still worthwhile to consider larger computational times for the 3-point correlators. Carrying 
out the complete low-mode averaging of the 2-point correlators takes about 35 minutes per 
configuration. The computation for the high modes with averaging over four time slices 
typically takes another 35 minutes, while averaging over all 48 time slices would typically 
take 7 hours. This means that choosing N ev = 100 and reducing the effect of the noisy high 
mode contribution is the optimal choice. 

6.1.2. Effect of low mode averaging for the 2-point correlator 

Although it is obvious that once we have the low eigenmodes we should carry out complete 
low-mode averaging, it would still be interesting to see how much gain in statistical accuracy 
we obtain with the low-mode averaging. For this purpose, we compare the results obtained 
with averaging over only four source time slices and those obtained with averaging over all 
time slices for the low-mode part. Figure [5] presents a comparison of the effective mass plots 
of the 2-point correlators. There, the low-mode contributions were averaged over four equally 
time slices separated and over 48 timeslices. There is obviously a significant improvement in 
the statistics. Considering the errors for t =5-10 we find that increasing N to from 4 to 48, 
the statistical error is further reduced by a factor of 2-3, which is effectively equivalent to 
having 4-9 more statistically independent configurations. This suggests that out of the 48 
time slices of 8-12 time slices are effectively statistically independent. 

6.1.3. Effect of to averaging for the high-mode part for the 2-point correlator 

We also studied the effect of t averaging for the high-mode contribution. In Fig. O we 
compare the data with N to = 4 and N to =48. It is seem that we obtain a clearer plateau 
for the data with Nt =48. In fact, the errors for the data with Nt = 4 are larger by a 
factor of approximately 1.5-2 than those with N to = 48 due to the reduction in the error 
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Fig. 5. Plots of the effective mass of the 2-point correlators for k = 0.1335 with low mode averaging 
over 4 time slices (top panel) and 48 time slices (bottom panel), with N r = 1 and N cv = 100. 



for the high mode part, which is not a negligible contribution. This is because, even though 
the high-mode contribution itself is small, its error is roughly equal to that of the low- mode 
parts. 

The data with N to = 48 in Fig. also exhibits a clean plateau which starts at t = 7, or 
perhaps even a smaller value. From this result, in the study of 3-point correlators reported 
below, we chose £^ = 7,8, and 9. 

6.1.4. N r dependence of the 2-point correlator 

It is also possible to reduce the error by increasing N r . However, as pointed out by 
the TrinLat group, increasing N r only reduces the error from the noisy estimator in the 
high mode, while increasing N to reduces the error from both the noisy estimator and the 
fluctuation of the gauge field. Since the numerical cost of computing the 2-point correlator 
is relatively small, we decided to use N to = 48 and N r = 1, which already gives sufficiently 
good results. We postpone the actual study of the N r dependence to a later subsection, 
where we study the 3-point correlators. 
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Fig. 6. Plots of the effective mass of the 2-point correlators for k = 0.1335 with complete low 
mode averaging. The high-mode part is averaged over 4 time slices (top) and 48 time slices 
(bottom), with N r = 1 and iV cv = 100. 



6.1.5. Low mode dominance in the 3-point correlator 

Now we turn to the 3-point correlators, for which we study the contribution of low-low, 
high-low, low-high and high-high modes. Figure [7] shows the contribution of each mode to 
the 3-point correlator for k = 0.1335, iV cv = 100, (N r ,N to ) = (1,4) and t A = 8. We find 
that for t < 9 (t + t A < 17) high-low and low-low are the two dominant contributions, with 
high- low larger for t < 4 (t + t A < 12) and low-low larger for t > 5 (t + t A > 13). 

While it is true that low-mode averaging does help to reduce the error, it is also important 
to reduce the error for the high-low contribution by carrying out more averaging for the high- 
modes. This is realized in two independent ways. One is to take a larger iV to for averaging 
the high- mode over the timeslices of the source points. The other is to increase N r . In the 
next subsections we examine the effects of both. 
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Fig. 7. Contributions of the low-low, high-low, low-high and high-high parts to the total 3-point 
function. 

6.1.6. (N r ,N tQ ) dependences of the 3-point correlator 

Figure [8] plots the (N r ,N to ) dependences of the 3-point correlator with t& = 8. We find 
that increasing N r from 1 to 4 reduces the statistical errors by a factor of approximately 
0.7, while increasing N to from 1 to 4 reduces the statistical errors for all t by a factor of 
approximately 0.5. Therefore, we find that for given computational resources, increasing N to 
is more efficient than increasing N r . This is as expected, since increasing N to reduces both 
the error from the gauge fluctuation and the error from the noisy estimator, while increasing 
N r reduces only the latter. 

6.1.7. tA dependence 

Figure M plots the time dependence of the 3-point /2-point ratio R(t, t^) for = 7, 8 and 
9. The results for t a = 7, 8 and 9 are consistent. The result of the data fit for tA = 9 is 
relatively poor and has large statistical error, while those for tA = 7 and 8 do not differ in 
precision. From this study, we conclude that tA = 8 is a reasonable choice for the productive 
run. We also find that the best fitting range of R(t, tA) is t + tA = 12 — 16(t = 4 — 8) . 

§7. Results 

In this section, we present the results from our productive run for 32 cofingurations on 
16 3 x 48 lattices with (3 = 6.0 using the 0(a)-improved Wilson fermion with csw = 1-769. 
We chose three values of the hopping parameters, k = 0.1335, 0.1340, and 0.1342. For the 
all-to-all propagators, our choices of the parameters were 

N ev = 100, N t0 = 48, N r = 1 (7-1) 
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Fig. 9. tA dependence of the 3-point/2-point correlator ratio for k = 0.1335 with (N r , N to ) = (1, 1) 
and N ev = 100. 
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Fig. 10. The ratio Rns f° r the 2-point correlation function with 32 configurations for k = 0.1335. 



7.1. 2-point correlators 

Figure fTUI displays the noise-to-signal ratio Rns = AC2(t) / C2{t) for the 2-point function. 
As shown in this figure the noise-to-signal ratio remains roughly in the range 1-3 % for the 
time range 1-2 fm. In this subsection we present our results for the 2-point correlators. 
Figure fTTl presents effective mass plots of the 2-point functions for k = 0.1335, 0.1340, and 
0.1342. Fitting the 2-point function with the fit range t = [7,11], we obtain the binding 
energy of the heavy-light meson with 0.2 % accuracy as shown in Fig. [121 The square in 
Fig. [12] represents the result of the Alpha collaboration® for the same value j3 = 6.0 on a 
16 3 x 32 lattice with the Schrodinger functional boundary conditions using on the order of a 
few hundred configurations. Although a simple comparison may not be meaningful, because 
of the different the kinematical setups used (for example, the volume, boundary conditions 
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Fig. 11. Plots of the effective mass of the 2-point function for k = 0.1335, 01340, and 0.1342, with 
(N r ,N to ) = (1,48). We observe a clear plateau for t £ [7, 11] 
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Fig. 12. Binding energies from 2-point correlators. The solid circles represent our data for k = 
0.1335, 0.1340, and 0.1342 with 32 configurations. The square represents the result from the 
Alpha collaboration with a few hundred configurations. 

and operators), it is nevertheless noteworthy that the all-to-all propagators yields much 
better precision with a much smaller number of configurations, iV con f = 32. Extrapolating 
linearly in (am n ) 2 , we obtain 

aE = 0.408 ± 0.002, (7-3) 

i.e., 0.5% accuracy in the chiral limit. 

7.2. 3-point correlators 

Figure [13] represents plots of the effective mass of the 3-point correlators. We find that 
the plateau is reached for t > 4. Figure [T41 displays the 3-point /2-point correlator ratios. 
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Fig. 13. Plots of the effective mass of the 3-point correlation function for re = 0.1335, 0.1340, and 
0.1342 with 4 to sources. 
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Fig. 14. 3-point/2-point correlator ratio R(t,t A ) for re = 0.1335, 0.1340, and 0.1342 with 4 t 
sources. We observe a clear plateau for t G [5, 9]. 



Applying constant fits of this ratio with fitting range t = [5, 9], we obtain 

R(t,t A ) = c^+f^ = 2-56(5) for « = 0.1335, 
= T^TT^TT = 2 - 55 ( 5 ) for K = °- 1340 ' 

= ^!f M \ = 2 - 5 °( 4 ) f ° r « = °- 1342 ' 

C 2 (t + EaJ 



(7-4) 
(7-5) 
(7-6) 



where each result has 2% accuracy. The bare lattice operator for the light-light axial vector 
current should be renormalized as 



Ai = 2kZ a (1 + b A m q )(Af + ca^P) 



(7-7) 
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Fig. 15. for k = 0.1335, 0.1340, and 0.1342, together with the result in the chiral limit obtained 
by linear extrapolation in (am ff ) 2 . 

The second term here does not contribute after summing over all spatial positions. The 
quantity Za has been determined nonperturbatively by the Alpha collaboration. For (3 = 6.0, 
they obtained Za = 0.79. The quantitly 6 a was also obtained by combining the results of by 
from the Alpha collaboration^ and bA — by from Battacharya et al. ® for f3 = 6.0. Using 
the above formula, the B*Bn coupling can be obtained as 

&o = 2kZ a (1 + b A m q )R. (7-8) 

Figure [15] displays our results the function of (am^) . Taking the chiral limit by 

extrapolating linearly in (am n ) 2 , we obtain g^ in the chiral limit as 

c/oo = 0.517 ±0.016, (for (3 = 6.0) (7-9) 

with 3% statistical error. This result can be compared with the result on 16 3 x 32 lattices for 
(3 = 6.0 using the 0(a)-improved Wilson fermion with csw = 1-769 by Abada et al.^ 1 Their 
result in the chiral limit is goo = 0.45(4), which is consistent with ours within 2a. We note 
that their results for finite quark masses have 4-6% statistical errors giving 9% accuracy in 
the chiral limit with on the order of 100 configurations, whereas ours have 2% statistical 
error for finite quark masses giving 3% accuracy in the chiral limit with 32 configurations. 
It is thus seen that are obtain a significant improvement in statistical accuracy using the 
all-to-all propagators. 

§8. Conclusions 

In this paper, we have studied the feasibility of the all-to-all propagator method for 
precise computation of the B*Bn coupling, g^. Determining 100 eigenmodes, we find that 
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Table II. Computational cost per configuration for 2-point and 3-point correlators for k = 0.1335. 









Z-pUlllL COIIclclLOI 
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ly t a 


low 


ni p*Vi 
lllgll 


total 
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4 




3D min 
'.j \j linn. 


.Xf) min 
'.j \j 111111 


20 


4 


7 min. 


35 min. 


40 min 


100 


1 


35 min. 


35 min. 


70 min 


100 


48 


35 min. 


7 h 25 min 


8 h 








3-point correlator 




iVcv 


Nt 


low-low 


low-high, high-low, high-high 


total 


100 


4 


40 min. 


~ 32 h 


~ 33 h 



the low-modes dominate the correlators, and therefore the low-mode averaging in the all-to- 
all propagator reduces the noise significantly. 

For reference we, list the numerical cost estimate using 1 CPU on our vector super- 
computer in Table [Til For 2-point heavy-light meson correlators in quenched QCD, 8 
hours are required per configuration using the all-to-all propagator with N cv = 100 and 
(N r ,Nt ) = (1,48), while this time becomes 10 hours if the time for the eigenmode calcu- 
lation is included. If the equivalent computational cost were used for the ordinary method 
with a single source timeslice, one would have to compute correlators with 64-80 times more 
configurations, which could reduce the error by a factor of 8-9. On the other hand, the 
errors of the 2-point correlator for 32 configuration are a factor of 0.08-0.10 smaller than 
those obtained with the ordinary method. Therefore, there is a marginal gain with the all- 
to-all propagator. We presume that the situation would be similar for 3-point correlators, 
although we have not made a comparison of the results obtained using all-to-all propagators 
and those obtained with the ordinary method. 

However, the essential aim of this study is not to show that we have a marginal gain 
with the all-to-all propagator in quenched QCD. What is important is that our result in 
quenched QCD suggests that there is a significant gain in unquenched QCD, in which case 
generating configurations is much more costly than the measurements. Also, the fact that 
only 32 independent gauge configurations are sufficient for obtaining very precise results with 
a few percent accuracy is encouraging. 
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In this appendix we give the definitions of the low-low, low-high, high-high, and high-low 
mode correlators, Cu, C\h, Chh, and Chi, where the first (second) index of C qiQ2 specifies 
the quark propagator from t A to t B * (from t B to t A )- The low-low contribution to the 3-pt 
functions is given by 



Appendix 



Cu(t B *,t A , t B ) 



V XT" 5ZS Tr {^Sl(x,t B *;t B )5(y - x)-f 5 r B * 



x [v®{x + r,t B *)®vW{z,t Ai ))r A 




(A-1) 



(A-2) 



where a = 1, 2, 3 and 




(A-3) 



z 




(A-4) 



X 



with r Aa = 7 a , r' B = r B 7 5 = -75, and r' B , = 75-T Bs = 757^. 
The low-high contribution to the 3-pt functions are given by 



Cih(t B *,tA,t B ) = — ^ y^Tr { S b (y, t B ; x, t B *)r B *Q (x + w, t B *; z, t A ) 




(A-5) 
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1 1 



i=l rj r,u> 

x [v®(x + r,t B *)®vW(z,t A )]r Aa 



X 



i i 



i=l r,j 



where 



(A-6) 

(A-7) 
(A-8) 



While we repeatedly use the symbols A^'(^) and B l J(t B *,t B ), we anticipate there will be no 
confusion, since they are valid only in the expression just before they are defined. 

In the case of Cm, we set the source of the noisy estimator to t B *. Since Q\ is hermitian, 
we have 



1 



QiO, z) = Q^z, x) f = — Wr]( Z ) ® V[ 



t _ 1 



This leads to 



C M (t B *,t A , t B ) = — Trjs'^y, t B ; f, tB*)r B *Qi(f + uT, 2, *a) 

x,y,z r,w 

x r Aa Q {z, t A ; y + w, t B )r B (f)(r)<f)(w) 
= y ^ X!^ E E E Tr {^( f > *b; *fl.W - 



(A-10) 



X 



J$(f+f,tB0®V$ T (3*A) 



(Ot, 



i i 



«=1 rj 



where 



B%(t B *,t B ) = ^uJ f)t (x,tfl)ffl 1 S' 6 (f,tB;tB.)rfl.77^(^tB. 



(A-ll) 

(A-12) 
(A-13) 
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Recall that _Tb* = j a , ^ = la, and t B = —1 (a = 1, 2, 3). 

The high-high part, C^ h is evaluated using the noisy estimator with a source vector put 
on the time slice t B : 

C hh (t B *,t A ,t B ) = — ^ ^ Tr ^S b (y, t B ; x, t B *)r B% Qi(x + w, t B * ;z,t A ) 

x,y,z r,w 

x r^Qx&tAiy + w,t B )f B (j)(r)(j)(w)} (A-14) 
1 1 



^Tr|s , 6 (y,t B ;f, t B *)r B *Qi{x + uj, t B *;z,t A ) 

r,i x,y,z r,w 

a [i>®(z,t A )®r,W(y + w,t B )] r B <j>{r)<t>{w). (A-15) 



V N r 

r,i x',y,z r,w 



x r A 

Using the source method, 

&[r] a (x> t B* 5 *a) = J] Ql (z, ^B* ; 3 tO-^^f*] *a) 

= 2 J D- 1 (^t B .;^t A )7 5 P 1 r j 4 a ^(z,t A ) (A-16) 

z 

can be obtained by solving the linear equation D\I/ = b, with the source vector b = 
j5Pir Aa ip$(z,t A ). Then C hh is obtained as 

c hh {t B *,t A ,t B ) = (A ' 1T) 

We must solve 3 x N tA linear equations per noise vector, where N tA is the number of sampled 
time slices t A . 
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